% built for five minute non-continous data.

%Items still to complete
% daily or nightly balancing
% Identifying storm events by sudden delta or change from mean to identify days when warping factors can't be calculated
% warping to better fit
% data tests and fixes
% Does it as written work for 15 minute data (No 10/17/22)
%TODO
% develop blanking routine from three sources
% Rain file
% Event File  ( see Prog_Dailypattern.m )
% flow values
% develop outlier routine


% 0.0 --- Initial Setup ---
clc
clearvars
clearvars -GLOBAL
format compact
warning('off')
tic
disp('********************** Starting New Run ***************************')


global PATH T NP DN PATTERN DEFINED Y M D H MN S 
RemoveEvents = 0;
Plotall = 1;


% disp(" ")
% disp("*********************************************************************************")
% disp(" John Barton 9/26/22")
% disp(" Currently identifies storm outliers.  Action items:")
% disp(" 1.) convert paseflow pattern to a daily flow pattern that balances before and after events.")
% disp(" 2.) use the mean value of before prior to the event, and the mean value of after  after the event")
% disp(" 3.) run the loop a second time with the ouliers (and missing values?) replaced with the above mean")
% disp(" 4.) sould sgolay be used outside of the outliers?")
% disp("*********************************************************************************")



% 1.0 *****------------ Load Path  -------------
disp('Load Data Directory')
[PATH] = uigetdir('*.csv', 'Select the Data Directory')
addpath (PATH)
cd (PATH)


% 1.1 ----------------  Load Flow Data  ---------------------------------------------
  % 1.1.1 ----------------  Load Flow File and interpret  ---------------------------------------------

  exist ("ModelOptFlowloaded.mat") 
  if exist ("ModelOptFlowloaded.mat") ~= 2
    [filename1, Path2] = uigetfile('*.csv;*.tsf','Select the FLOW DATA File')
    opts = detectImportOptions(filename1, 'FileType','text')
    [ fp1 , fn1 , ext1 ] = fileparts( filename1 )
    if ext1 == ".tsf"
        disp (" ")
        disp (strcat("There are ", num2str(opts.DataLines(1)-1), " header lines in the data file." ))
        fid = fopen(filename1);     
        frewind(fid)
        Varsin = 0;
        for kh0 = 1:opts.DataLines(1)-1
            line = fgetl(fid);
            disp(strcat("    ", line))
            FQt= strfind(line, "Date/Time");
            if FQt == 1
                Varsin = kh0;
                varNames =  line;
            end
        end
        if Varsin == 0
            Varsin = input ("What is the row number of the Header Line listed above which contains the Variable Names?   ")
            frewind(fid)
            for kh1 = 1:Varsin
                line = fgetl(fid);
                varNames =  line;
            end
        end
        disp(" ")
        disp("The Variable Names are ...")
        disp(strcat("     ", varNames))

        %Find the order of Flow, Velocity, and Level columns
        clear VarOrder
        %Flow
        FQt = strfind(varNames, "Edited Flow");
        if ~isempty(FQt)
            disp("Using Edited Flow column.")
        elseif isempty(FQt)
            FQt = strfind(varNames, "Flow");
        elseif isempty(FVt)
            FQt = strfind(varNames{kh}, "flow");
        end
        if length(FQt) > 1
            disp("WARNING: multiple Flow columns.  Using the first one.")
            FQt = min(FQt);
        end
        VarOrder.num{:,1} = FQt;
        VarOrder.name{:,1} = 'Flow';
        
        % Velocity
        FVt = strfind(varNames, "Edited Velocity");
        if ~isempty(FVt)
            disp("Using Edited Velocity column.")
        elseif isempty(FVt)
            FVt = strfind(varNames, "Velocity");
        elseif isempty(FVt)
            FVt = strfind(varNames{kh}, "velocity");
        end
        if length(FQt) > 1
            disp("'WARNING: multiple Velocity columns.  Using the first one.")
            FVt = min(FVt);
        end
        VarOrder.num{:,2} = FVt;
        VarOrder.name{:,2} = 'Velocity';
        
        FDt = strfind(varNames, "Edited Level");
        if ~isempty(FDt)
            disp("Using Edited Level column.")
        elseif isempty(FDt)
            FDt = strfind(varNames, "Level");
        elseif isempty(FDt)
            FDt = strfind(varNames{kh}, "Level");
        elseif isempty(FDt)
            FDt = strfind(varNames, "Edited Depth");
            if ~isempty(FDt)
                disp("Using Edited Depth column.")
            end
        elseif isempty(FDt)
            FDt = strfind(varNames, "Depth");
        elseif isempty(FDt)
            FDt = strfind(varNames{kh}, "Depth");
        end
        if length(FQt) > 1
            disp("'WARNING: multiple Level columns.  Using the first one.")
            FDt = min(FDt);
        end
        VarOrder.num{:,3} = FDt;
        VarOrder.name{:,3} = 'Level'

        disp(strcat("The variable name for the first column will be changed to Date."))
        [~, order] = sort(cell2mat(VarOrder.num));
        Orderednames = VarOrder.name(order);
        opts.VariableNames{1} = 'Date';
        for kord = 1:length(cell2mat(VarOrder.num))
            opts.VariableNames{kord+1} = Orderednames{kord};
        end
        opts = setvaropts(opts,'Date','DatetimeFormat','MM/dd/uuuu HH:mm');
        T = readtable(filename1,opts); % T.Date, T.Flow, T.Depth, T.Velocity

    elseif min(ext1 == '.csv')
        disp ("Because this is a .csv file, it is assumed to be a seven row header." )
        fid = fopen(filename1);
        opts.DataLines(1) = 8;
        disp(' ')
        disp (strcat("The first three rows of the file are " ))
        for kh = 1:3
            linetemp =fgetl(fid);
            disp(linetemp);
            if kh>1
                varNamesT{kh-1} =  split ((linetemp),",");
            end
        end

        for kt = 1:length(varNamesT{1})
            if kt == 1
                varNames{kt} = 'Date';
                disp(' ')
                disp(strcat("The variable name for the first column has been set to Date."))
            else
                varNames{kt} = strcat(varNamesT{1}{kt}, '_', varNamesT{2}{kt});
            end
            opts.VariableNames{kt} = varNames{kt};
        end
        disp (strcat("The variable names will be initially loaded to a table as:" ))
        varNames

        %    opts = setvaropts(opts,'Date','DatetimeFormat','M/d/uu H:mm')
        T = readtable(filename1,opts);
        T.Flow = T.edit_flow;
        T.Depth = T.edit_level;
        T.Velocity = T.edit_velocity;
    end



   % 1.1.2 ----------------- Add additional variable columns with zeros to the table  --------------------
    NP = length(T.Flow);
    WT = width(T);
    DN = datenum(T.Date);

    Outlier = zeros(NP,1); % New Column in Table
    T = addvars(T, Outlier); % T.Date, T.Flow, T.Depth, T.Velocity, T.Outlier
    clear Outlier
    RawBaseFlow = zeros(NP,1); % New Column in Table
    T = addvars(T, RawBaseFlow); % T.Date, T.Flow, T.Depth, T.Velocity, T.Outlier, T.BaseFlow
    clear RawBaseFlow
    BaseFlow = zeros(NP,1); % New Column in Table
    T = addvars(T, BaseFlow); % T.Date, T.Flow, T.Depth, T.Velocity, T.Outlier, T.BaseFlow
    clear BaseFlow
    StormFlow = zeros(NP,1); % New Column in Table
    T = addvars(T, StormFlow); % T.Date, T.Flow, T.Depth, T.Velocity, T.Outlier, T.BaseFlow, T.StormFlow
    clear StormFlow
    disp(' ')
    disp (strcat("Additional columns have been added to the table:" ))
    varNamesUpdated = T.Properties.VariableNames

    % 1.1.3 ---------------- Some Data Corrections  --------------------
    % Replace 999
    T.Flow(find(T.Flow==-999))=0;
    if max(strcmp('Depth',T.Properties.VariableNames))
        T.Depth(find(T.Depth==-999))=0;
    end
    if max(strcmp('Velocity',T.Properties.VariableNames))
        T.Velocity(find(T.Velocity==-999))=0;
    end

    % Determine time step of loaded data
    Timesteps = (DN(2:end)-DN(1:end-1));
    TSperDay = round(1/median(Timesteps));
    Timestep = 1/TSperDay;

    if (max(Timesteps) - min(Timesteps)) < 0.00000001
        disp(strcat("All Timesteps equal ",  num2str(Timestep*24*60), " minutes."))
    else
        disp( "WARNING: There are probably Missing or Duplicate Time Stamps.  Missing will be added and filled with zeros.  Duplicates will be deleted.")
        format long
        MedianTimestep = median(Timesteps);
        MaxTimestep = max(Timesteps);
        MinTimestep = min(Timesteps);
        format short
        if MaxTimestep ~= Timestep || MinTimestep ~= Timestep
            disp(strcat("Maximum Timestep ",  num2str(MaxTimestep*24*60), " minutes."))
            disp(strcat("Minimum Timestep ",  num2str(MinTimestep*24*60), " minutes."))
        end

        % Add missing datestamps
        for k = 2:NP
            if round(TSperDay*(DN(k)-DN(k-1)),  4) ~=1
                Num_pointsadd = round(-1+TSperDay*(DN(k)-DN(k-1)))
                if Num_pointsadd == -1  %delete duplicates
                    disp( "Duplicate ")
                    disp(string(T.Date(k-1,:)))
                    disp(string(T.Date(k,:)))
                    T = [T(1:k-1,:); T(k+1:end,:)];
                    disp ('  Duplicate time stamp deleted.')
                else
                    disp(strcat( "Prior to missing data ", string(T.Date(k-1,:)), "   After missing data ", string(T.Date(k,:))))
                    T=T([1:k-1+Num_pointsadd,k:end],:);
                    for k3 = 1:Num_pointsadd
                        ra = k+k3-1;
                        T(ra,2:width(T)) = num2cell(0);
                        T(ra,1) = cellstr(datestr(DN(k-1)+k3*Timestep,'mm/dd/yyyy HH:MM' ));
                        T.Outlier(ra) = 1;
                    end
                end
            end
        end
        NP = length(T.Flow);
        DN = datenum(T.Date);
        Timestep = ((DN(end)-DN(1))/(length(DN)-1));
        disp( " ")
        disp("******* End of Duplicates and Missing *******")
    end
    TSperDay = 1/Timestep;

    %Remove leading rows with no data just a timestamp

    %Remove ending rows with no data just a timestamp


  %   %1.1.4 ---------------  Save data for quick load during debugging -----------
  %   disp( " ")
  %   save('ModelOptFlowloaded.mat', 'T', 'PATH', 'filename1', 'Path2', '-mat');
  %   disp(strcat("Saved file   ", PATH, "ModelOptFlowloaded.mat"))
  end
  % 
  %     % %clear everything and load the file. It will skip to this if the ModelOptFlowloaded.mat file exists.
  %     clearvars
  %     global PATH T NP DN PATTERN DEFINED Y M D H MN S
  %     RemoveEvents = 0;
  %     Plotall = 1;
  %     pwd
  %     load ('ModelOptFlowloaded.mat', '-mat'); % loads T, PATH, filename1, Path2
  %     NP = length(T.Flow);
  %     DN = datenum(T.Date);
  %     Timestep = ((DN(end)-DN(1))/(length(DN)-1));
  %     TSperDay = 1/Timestep;


  %1.1.5 ---------------  OPTIONAL: For Debugging: reduce the size of the data set -------------------
    % lengthofdatadesired = 50000;
    % if length(T.Flow) > lengthofdatadesired
    %     format long
    %     Cutoff = floor((length(T.Flow)-lengthofdatadesired)/2 )
    %     format short
    %     T([1:Cutoff],:) = []; T([end-Cutoff:end],:) = [];
    %     disp(" "), disp(strcat("The length of the data set has been truncated to ", num2str(length(T.Flow)), " points for code development."))
    %     NP = length(T.Flow)
    %     DN = datenum(T.Date);
    % end




  % 1.1.6 ----------------  Display Details from the Flow Data ----------------------
    [Y,M,D,H,MN,S] = datevec(DN);
    [Day1num, Day1nam] = weekday(DN(1));
    [Day1endnum, Dayendnam] = weekday(DN(end));
    disp(" ")
    disp(strcat(['First time stamp is ' Day1nam(1,1:3) ' ' datestr(T.Date(1), 'mmm dd, yyyy HH:MM')]))  % for debug
    disp(strcat(['Last time stamp is ' Dayendnam(1,1:3) ' ' datestr(T.Date(end), 'mmm dd, yyyy HH:MM')]))  % for debug
    disp(strcat("The data period spans ", num2str(round(DN(end)-DN(1),4)) , " days."))   % for debug
    disp(strcat("The number of data points is equivalent to ", num2str(round((Timestep*NP),4)) , " days."))   % for debug
    disp(strcat("There appear to be ", num2str(round(1+ (DN(end)-DN(1))/Timestep - NP), 2), " missing time stamps."))   % for debugging
    FirstMondp = find( (weekday(DN)==2) .* (H==0) .* (MN<5), 1 );




%1.2 ----------------  Load Rain Data ----------------------
    %1.4 ---------------  Save Rain data for quick load during debugging -----------

    if exist ("ModelOptRainloaded.mat") == 2
        pwd
        load ('ModelOptRainloaded.mat', '-mat') % loads T
  
    else
        disp( "Loading Rain")
        [Rainname3, Path3] = uigetfile('*.dat;*.tsf', 'Select the RAIN DATA File')
        optsR = detectImportOptions(Rainname3, 'FileType','text')
        [ fp1 , fn1 , ext1 ] = fileparts( Rainname3  )

        if ext1 == ".tsf"

            %%%
            disp (" ")
            disp (strcat("There are ", num2str(optsR.DataLines(1)-1), " header lines in the Rain file." ))
            fid = fopen(Rainname3);
            frewind(fid)
            VarsinR = 0;
            for kh0 = 1:optsR.DataLines(1)-1
                line = fgetl(fid);
                disp(strcat("    ", line))
                FRt = strfind(line, "Date/Time");
                if isempty(FRt)
                    FRt = strfind(line, "Date");
                elseif isempty(FRt)
                    FRt = strfind(line, "date");
                end
                if FRt == 1
                    VarsinR = kh0;
                    varNamesR =  line;
                end
            end
            if VarsinR == 0
                VarsinR = input ("What is the row number of the Header Line listed above which contains the Variable Names?   ")
                frewind(fid)
                for kh1 = 1:VarsinR
                    line = fgetl(fid);
                    varNamesR =  line;
                end
            end


            disp(" ")
            disp (strcat("The variable names in the file appear to be:" ))
            disp(strcat("     ", varNamesR))


            disp(strcat("The variable name for the first column will be changed to Date and the second to Rain."))
            optsR.VariableNames{1} = 'Date';
            optsR.VariableNames{2} = 'Rain';

            %optsR = setvaropts(opts,'Date','DatetimeFormat','MM/dd/uuuu HH:mm');
            R = readtable(Rainname3,optsR); % R.Date, R.Rain
            if max(R.Date(2:end) - R.Date(1:end-1)) ~= min(R.Date(2:end) - R.Date(1:end-1))
                disp("There appears to be a time step issue.  It might be missing timesteps or")
                disp("Matlab might have misinterpreted the date format. ")
                disp("This would require modifications to the code.")
                disp("Below is some information to help figure it out.")

                Max_Rain_Timestep = max(R.Date(2:end) - R.Date(1:end-1))
                Min_Rain_Timestep  = min(R.Date(2:end) - R.Date(1:end-1))
                dtemp = R.Date   % datetime array

                dtemp.Format
            end

        elseif min(ext1 == '.dat')
            disp ("Because this is a .dat file, it is assumed to have no header." )
            fid = fopen(Rainname3);
            disp(' ')
            disp (strcat("The first three rows of the file are " ))
           
                    optsR.VariableNames{1} = 'RainGauge';
                    optsR.VariableNames{2} = 'Ryear';
                    optsR.VariableNames{3} = 'Rmon';
                    optsR.VariableNames{4} = 'Rday';
                    optsR.VariableNames{5} = 'Rhour';
                    optsR.VariableNames{6} = 'Rmin';
                    optsR.VariableNames{7} = 'R_Depth';


            disp (strcat("The variable names will be initially loaded to a table as:" ))
            optsR.VariableNames

            %    opts = setvaropts(opts,'Date','DatetimeFormat','M/d/uu H:mm')
            R = readtable(Rainname3,optsR);
            %keyboard
        end
    end

    









% 2.0 *****----------------  Daily patterns ----------------------
% 2.1 ---------------  hold out middle zeros --------------- 
    Zeroindices = find (T.Flow ==0);

% 2.2 ---------------  Basic Daily Patterns --------------- 
    % Smoothing
    Sm = 48; % smoothing factor.  Use 24 to look two hours each way in 5 min data
    Qs = smoothdata (T.Flow, 'sgolay', Sm, 'degree', 2);  % smooth data before making a daily pattern
    plot (Qs)
    % Daily Patterns
    whos global
    whos Qs DN
     F_DailyPatternGenerate_1day (Qs, DN)
     F_DailyPatternGenerate_2day (Qs, DN)
     F_DailyPatternGenerate_diff (Qs, DN)
     F_DailyPatternGenerate_Median (Qs, DN)
     for Type = 1:4  % 4 types of Daily Patterns
         Dailypatterns{Type} = F_baseTScreate(DN,Type);
     end % for Type = 1:4
    
     %Select the daily Patterns to use based on the length of the data series
     if DN(end)-DN(1) < 22  %use one day if less than 22 days of data
         SPat = 4;
     elseif  DN(end)-DN(1) < 36  %use weekday weekend if less than 36 days of data
         SPat = 3;
     else  % with 36 or more days of data use the 7-day by median
         SPat = 2;
     end
%SPat = 4

     disp(" ")
     disp(strcat("Selected Pattern is ",  PATTERN(SPat).LongName))
         disp(" ")

 T.RawBaseFlow= Dailypatterns{SPat};


 % coffset calc

 % disp(strcat("For testing the 29.5 in, diameter of MC-EB-120 is being used. "))
 % Diameter = 29.5/12; % feet
 % Cslope = offset_equation(c_offset, V, 0.2/12, Dia, depth)
 % c_offset = calculate_c_offset(T.Velocity, Cslope, Diameter, T.Depth)

 
 
 
 
 %3.0-Alt Load Rainfall





% 3.0 *****----------------  Outliers and Storm Events ----------------------
% 3.1 ----------------   Outliers  ---------------- 
    % Outliers
    Window1 = 288; % half a day each way
    [DPFit_byPoint] = DPatternFit(Qs, Dailypatterns{SPat},Window1) ; % half a day each way);

    afterOutlier = 12*24;  %how much to add to outliers after the last found point
    beforeOutlier = 12*4; %how much to add to outliers before the firstlast found point
    FillbetweenOutliers = 12*12; % half a day each way
    [OutliertoRemove, ~] = OutlierRules(DPFit_byPoint{4}, 0,0,0); 

    %storm outliers
    afterOutlier = 12*24;  %how much to add to outliers after the last found point
    beforeOutlier = 12*4; %how much to add to outliers before the firstlast found point
    FillbetweenOutliers = 12*12; % half a day each way
    [StormOutlier_Ind, NonStormOutlier_Ind] = OutlierRules(DPFit_byPoint{4}, afterOutlier, beforeOutlier, FillbetweenOutliers); 
    DN_NonStormOutlier = DN(NonStormOutlier_Ind);
    DN_StormOutlier = DN(StormOutlier_Ind);
    Outlier_Val = DPFit_byPoint{4}(StormOutlier_Ind);
    
%     Outlier_Flow = T.Flow(Outlier_Ind);
%     NonOutlier_Flow = T.Flow(NonOutlier_Ind);

    % Determine balance points before storm
    % AveBeforeStorm contins information about the identified period just before each outlier or outlier group.
    [AveBeforeStorm]  = f_avebeforestorm(Qs,  Dailypatterns{SPat});


% 2.4 --- BaseFlow ---
    CorrectionToBaseflow = interp1( [DN(1); DN(AveBeforeStorm.meantimeindex); DN(end)],...
        [AveBeforeStorm.BaseflowCorrection(1) AveBeforeStorm.BaseflowCorrection AveBeforeStorm.BaseflowCorrection(end)],...
        DN,'linear');
    BaseFlowTemp = Dailypatterns{SPat} + CorrectionToBaseflow;

    %Baseflow Options
    BFopt = 1;
    if BFopt == 1
        T.BaseFlow =  BaseFlowTemp;
    elseif BFopt ==2 % this option will use sgolay outside of the storms.
        T.BaseFlow = (Qs .* (T.Outlier==0)) + (BaseFlowTemp .* T.Outlier);
    end


% 2.5 --- Correct Flow for Missing data
    MissingData = find(T.Flow ==0); % these are points with time Stamps.
    if length(MissingData) > 0
        disp (' ')
        disp('Not Yet Written code for missing data with datestamps must ')
        disp('   1.) Find the points on either side, ')
        disp('   2.) Correct the baseflow to have the proper slope, ')
        disp('   3.) Replace the flow with corrected slope baseflow.')
        disp('   4.) Backcalculate the missing data channel (level or velocity) ')
        disp ('In the meantime, it just uses the baseflow')

        T.Flow(MissingData) = T.BaseFlow(MissingData);
        DNreplaced = DN(MissingData);
        Flowreplaced = T.Flow(MissingData);
    end

% 2.6 --- 
    T.StormFlow = T.Flow-T.BaseFlow;
    [~, DailyStormavebypoint, ~] = DailyAverage(T.StormFlow);
    [DailyAve, Dailyavebypoint, MoveAve24cen] = DailyAverage(Qs);
    Mindays = find(islocalmin(DailyAve),1)


%*********************************************************************************************

if Plotall == 1
    figure(1)
    disp('Plotting Figure 1')
    clf
    title("Raw Flow Data and Smoothed with Savitsky-Golay")

    hold on
    plot(DN, T.Flow)
    plot(DN,Qs)
    legend (' Raw Data', 'S-G Smoothed data')
    hold off
    % export current axes to JPEG at 300 DPI
    graphname = 'Figure 1 - Raw Flow Data and Smoothed with Savitsky-Golay.jpg';
    ax = gca;
    exportgraphics(ax, graphname, "Resolution", 300);



    figure (2)
    disp('Plotting Figure 2')
    clf
    title(strcat("Raw Flow Data and Daily Pattern from ", PATTERN(SPat).name))
    hold on
    plot (DN,T.Flow)
    plot (DN, Dailypatterns{SPat})
    legend ('Raw Data', 'From Selected Daily Pattern - not balanced')
    hold off
    graphname =  strcat("Figure 2 - Raw Flow Data and Daily Pattern from ", PATTERN(SPat).name,".jpg");
    ax = gca;
    exportgraphics(ax, graphname, "Resolution", 300);


    f = figure(3)
    disp('Plotting Figure 3')
    clf(f)
    title("The four different calculations of Daily Pattern - In Parallel.")
    hold on
    for Type = 1:4  % 4 types of Daily Patterns
        ax(Type) = subplot(4,1,Type)
        plot(PATTERN(Type).as7day)
        title(PATTERN(Type).LongName)
    end % for Type = 1:4
    hold off
    graphname = 'Figure 3 - The four different calculations of Daily Pattern - In Parallel.jpg';
   drawnow
    exportgraphics(f, graphname, "Resolution", 300);


    figure (4)
        disp('Plotting Figure 4')
    clf
    title("The four different calculations of Daily Pattern.")
    hold on
    for kty = 1:4
        %     plot (DN,Dailypatterns{kty})
        plot (Dailypatterns{kty}(1:12*24*7+1))
    end
    legend ('7-day Difference', '7-day Median', ' Weekday-weekend 6:00 PM', '1-day Median')
    hold off
    drawnow
    graphname = 'Figure 4 - The four different calculations of Daily Pattern - Overlaid.jpg';
    ax = gca;
    exportgraphics(ax, graphname, "Resolution", 300);


    % figure (5)
    % disp('Plotting Figure 5')
    % clf
    % title("Daily Average by Point and Daily Pattern Average.")
    % hold on
    % plot (DN,Dailyavebypoint)
    % plot (DN,Dailyavebypoint_Pat)
    % legend ('Daily Average of Raw Data', 'Daily Average of Daily Pattern')
    % hold off

    figure (6)
    disp('Plotting Figure 6')
    clf
    hold on
    plot (DN,DPFit_byPoint{4}')
    plot (DN(find(T.Outlier)), DPFit_byPoint{4}(find(T.Outlier)), 'r.')
    hold off
    graphname = 'Figure 6 - Outliers.jpg';
    ax = gca;
    exportgraphics(ax, graphname, "Resolution", 300);


    figure (7)
    disp('Plotting Figure 7')
    clf
    title("Outliers identifed by daily pattern matching on both sides, but not in the middle.")
    hold on
    plot (DN,T.Flow)
    plot (DN(find(T.Outlier)), T.Flow(find(T.Outlier)), 'r.')
    legend ('Raw Flow Data', 'Outliers')
    hold off
    graphname = 'Figure 7 - Outliers identifed by daily pattern matching on both sides but not in the middl.jpg';
    ax = gca;
    exportgraphics(ax, graphname, "Resolution", 300);

    figure(8)
    disp('Plotting Figure 8')
    clf
    title("Flow by outlier or not")
    hold on
    plot(DN, T.Flow)
    plot(DN(find(T.Outlier)), T.Flow(find(T.Outlier)), 'r.')
    plot(DN(find(T.Outlier==0)), T.Flow(find(T.Outlier==0)), 'g.')
    legend ('Raw Data', 'Outlier Flow', 'NonOutlier Flow')
    hold off
    graphname = 'Figure 8 - Flow by outlier or not.jpg';
    ax = gca;
    exportgraphics(ax, graphname, "Resolution", 300);


    figure(9)
    disp('Plotting Figure 9')
    clf
    hold on
    plot (DN, T.Flow)
    plot (DN, CorrectionToBaseflow)
    legend (' Raw Data', 'NDPfit Adj')
    hold off
    graphname = 'Figure 9 - Raw Flow Data and Baseflow Adjustment.jpg';
    ax = gca;
    exportgraphics(ax, graphname, "Resolution", 300);


    figure (10)
    disp('Plotting Figure 10')
    clf
    title("Raw Flow with Savitsky-Golay and Balanced Daily Flow (currently selected days)")
    hold on
    plot (DN,T.Flow,'LineWidth',2.0)
    plot(DN,Qs,'LineWidth',2.0)
    plot (DN, Dailypatterns{SPat})
    plot(DN, CorrectionToBaseflow)
    plot(DN(AveBeforeStorm.meantimeindex),AveBeforeStorm.BaseflowCorrection, 'ro')
    legend ('Raw Data', 'SGolay', 'Baseflow Unadjusted', 'Baseflow Adjustment')
    hold off
    graphname = 'Figure 10 - Raw Flow Data and Smoothed with Savitsky-Golay (currently selected days).jpg';
    ax = gca;
    exportgraphics(ax, graphname, "Resolution", 300);



    % figure (11)
    % disp('Plotting Figure 11')
    % clf
    % title("Baseflow Adjustment")
    % hold on
    % plot (DN, Dailypatterns{SPat})
    % legend ('Baseflow Unadjusted')
    % hold off


    figure (12)
    disp('Plotting Figure 12')
    clf
    title("BaseflowTemp")
    hold on
    plot (DN, T.Flow)
    plot (DN, BaseFlowTemp)
    legend ('BFlow', 'BaseFlowTemp')
    hold off
    graphname = 'Figure 12 - BaseflowTemp.jpg';
    ax = gca;
    exportgraphics(ax, graphname, "Resolution", 300);


    figure (13)
    disp('Plotting Figure 13')
    clf
    title("Baseflow")
    hold on
    plot (DN, T.Flow, 'b')%, 'LineWidth',1.1)
    plot (DN, T.BaseFlow, 'r','LineWidth',1.2)
    % plot(DNreplaced, Flowreplaced, 'r.')
    xlmin = floor(DN(2000));
    xl = [xlmin xlmin+20];
    xlim(xl);
    %ylim([ 0 2.5])
    legend ('Flow', 'T.BaseFlow')
    hold off
    graphname = 'Figure 13 - Baseflow.jpg';
    ax = gca;
    exportgraphics(ax, graphname, "Resolution", 300);


    figure (14)
    disp('Plotting Figure 14')
    clf
    title("Storm Flow")
    hold on
    plot (DN, T.StormFlow)
    legend ('Storm Flow')
    hold off
    graphname = 'Figure 14 - Storm Flow.jpg';
    ax = gca;
    exportgraphics(ax, graphname, "Resolution", 300);


    figure (15)
    disp('Plotting Figure 15')
    clf
    title("Daily Mean Storm Flow")
    hold on
    plot (DN, DailyStormavebypoint)
    legend ('Daily Mean Storm Flow')
    hold off
    graphname = 'Figure 15 - Daily Mean Storm Flow.jpg';
    ax = gca;
    exportgraphics(ax, graphname, "Resolution", 300);


    if  sum(ismember(T.Properties.VariableNames,'Depth'))  && sum(ismember(T.Properties.VariableNames,'Velocity'))
        figure (16)
        disp('Plotting Figure 16')
        clf
        title("Level-Velocity")
        hold on
        yyaxis left
        plot (DN, T.Velocity, 'g')%, 'LineWidth',1.2)
        yyaxis right
        plot (DN, T.Depth, 'k')%, 'LineWidth',1.1)
        legend ('Level', 'Velocity')
        xlim(xl)
        hold off

        figure (17)
        disp('Plotting Figure 17')
        clf
        hold on
        title("Scattergraph")
        plot (T.Depth, T.Velocity, 'b.')
        plot (T.Depth(find(DN>xl(1)&DN<xl(2))), T.Velocity(find(DN>xl(1)&DN<xl(2))), 'r.')
        xlabel('Level')
        ylabel('Velocity')
        hold off
    else
        figure (16); close
        figure (17); close
    end

    disp(" ")
    toc
    disp(" ")
    disp("** End of Program **")

end


%for Plotting

% for k = 15000:2:length(Qs)
%     imin = max([1 k-Window1]);
%     imax = min([k+Window1-1 length(Qs)]);
%     cen = (imax-imin+1)/2;
%     figure (20)
%     clf
%     title("Flow and Daily Pattern normalized in window ")
%     hold on
%     plot (Qs(imin:imax)/mean (Qs(imin:imax)),'g')
%     plot  (Dailypatterns{SPat}(imin:imax)/ mean (Dailypatterns{SPat}(imin:imax)), 'b')
%     plot ([cen cen], [0 max(Qs(imin:imax)/mean (Qs(imin:imax)))], 'r-')
%     hold off
% 
%     figure (21)
%     clf
%     title("slope of dp subtracted residuals")
%     hold on
%     plot (DPFit_byPoint{1}(imin:imax), 'm-')
%     plot (DPFit_byPoint{3}(imin:imax), 'c-')
%     plot ([cen, cen],[0, max(DPFit_byPoint{1}(imin:imax))], 'r-')
% 
%     hold off
%     pause(.01)
% end

%*************************************************************************
function  balancedaily
end


%*************************************************************************
function balanceweekly
end

%*************************************************************************
function [DailyAve, Dailyavebypoint, MoveAve24cen] = DailyAverage(TS)
   global NP DN
   %TS could be any time series of level, velocity, or flow
   Dayofseries = fix(DN);
   Days = (min(Dayofseries):max(Dayofseries))';
   DailyAve = zeros(length(Days),1);
   Dailyavebypoint =  zeros(NP,1);
   for k9 = 1:length(Days)
       DailyAve(k9,1) = mean(TS((Dayofseries)==Days(k9)));
       Dailyavebypoint( fix(DN)==Days(k9),1) = DailyAve(k9);
   end
   Window1 = 144; % half a day each way
   MoveAve24cen = zeros(size(TS));
   for k=1:NP
       indexmin = max([1 k-Window1]);
       indexmax = min([k+Window1-1 NP]);
       MoveAve24cen(k) = mean(TS(indexmin:indexmax));
   end
end

%*************************************************************************
function [DPFit_byPoint] = DPatternFit(TS_Smooth, cur_dp, Window)

    for k=1:length(TS_Smooth)
        imin = max([1 k-Window]);
        imax = min([k+Window-1 length(TS_Smooth)]);
    
        % Forward
        TS_m = TS_Smooth(k:imax) /mean (TS_Smooth(imin:k));
        cdpm = cur_dp(k:imax) / mean (cur_dp(imin:k));
        DPFit_byPoint{1}(k) = (sum((TS_m-cdpm).^2))^0.5; %forward
    
        %centered
        TS_m = TS_Smooth(imin:imax) /mean (TS_Smooth(imin:k));
        cdpm = cur_dp(imin:imax) / mean (cur_dp(imin:k));
        DPFit_byPoint{2}(k) = (sum((TS_m-cdpm).^2))^0.5; %centered
    
        %backward
        TS_m = TS_Smooth(imin:k) /mean (TS_Smooth(imin:k));
        cdpm = cur_dp(imin:k) / mean (cur_dp(imin:k));
        DPFit_byPoint{3}(k) = (sum((TS_m-cdpm).^2))^0.5; %backward
    
        %Min
        DPFit_byPoint{4}(k) = min ([DPFit_byPoint{1}(k) DPFit_byPoint{2}(k) DPFit_byPoint{3}(k)]);


        


    end 
end


%*************************************************************************
function [Outlier_Ind, NonOutlier_Ind] = OutlierRules(fitbypoint,afterTS, beforeTS,FillWindow) 
    global T NP DN
    %Identify Outliers

%add time before and after the identified outliers

    
    OUTLIERStemp = isoutlier(fitbypoint);

    for kout = afterTS+1:NP-beforeTS
        if max(OUTLIERStemp(kout-afterTS:kout+beforeTS)) > 0
            T.Outlier(kout) = 1;
        end
    end


    % make outliers of small sections of data between other outliers
    
    OUTLIERStemp = T.Outlier;
    for k1=2:NP-1
        imin = max([1 k1-FillWindow]);
        imax = min([k1+FillWindow-1 NP]);
        if  T.Outlier(k1-1) ~=1 && T.Outlier(k1) ~=1 && sum(OUTLIERStemp(imin:k1-2))>0 && sum(OUTLIERStemp(k1+1:imax))>0 
            T.Outlier(imin:imax) = 1;
        end
    end
    if sum(T.Outlier(1:FillWindow)) > 0
        T.Outlier(1:FillWindow) = 1;
    end

    Outlier_Ind = find(T.Outlier);
    NonOutlier_Ind = find(T.Outlier==0);


end

%*************************************************************************
function [AveBeforeStorm]  = f_avebeforestorm( TS_Smooth,  cur_dp)
    global T NP DN
   % DPFit_byPoint_Adj = zeros(NP,1);
    count = 0;
    Window1 = 144; % half a day each way
    for k1=Window1:NP-1
        imin = max([1 k1-Window1]);
        imax = min([k1+Window1-1 NP-1]);
        if T.Outlier(imax+1) >0 && sum(T.Outlier(imin:imax)) == 0  % identifies the last 24 hours prior to a storm
            % finds the data point of the center of the last 24 hr period prior to the outlier
            TS_m = TS_Smooth(imin:imax) /mean(TS_Smooth(imin:imax));
            cdpm = cur_dp(imin:imax) / mean(cur_dp(imin:imax));
           % DPFit_byPoint_Adj(imin:imax) = cur_dp (imin:imax) - (-mean(TS_Smooth(imin:imax)) + mean (cur_dp(imin:imax)));
    
            count = count +1;
            ABS_tempcount(count) = count;
            ABS_tempmeantimeindex(count) = k1; % index of the time of the center point of the day (24-hr period) before the storm (outlier)
            ABS_tempmeantime(count) = DN(k1); % time of the center point of the day (24-hr period) before the storm (outlier)
            ABS_temptimeindex(count) = imax; % index of the time of the last data point  before the storm (outlier)
            ABS_temptime(count) = DN(imax); % time of the last data point  before the storm (outlier)
            ABS_tempTSmean(count) =mean(TS_Smooth(imin:imax)); % mean of the Flow or time series prior to the storm
            ABS_tempcdpmean(count) =mean(cur_dp(imin:imax)); % mean of the calculated daily pattern prior to the storm
            ABS_tempBaseflowCorrection(count) = ABS_tempTSmean(count) - ABS_tempcdpmean(count) ; % TS_Mean-cdpmean

            if count == 1
                GWIslope(count) = 0; %slope of the GWI increase or decrease
            else
                GWIslope(count) = (ABS_tempcdpmean(count) - ABS_tempcdpmean(count-1))...
                    /  (ABS_tempmeantime(count) - ABS_tempmeantime(count-1));
            end
        end
    end
    disp(strcat("Initial scan identified ", num2str(count), " storms starts."))

    % This section will delete the prestorm (preoutlier groups) points if it shows a sudden increase in GWI.
    GWI_IncreaseOK = find(GWIslope < mean(abs(GWIslope)));
    ABS_tempcount = ABS_tempcount(GWI_IncreaseOK);
    ABS_tempmeantimeindex = ABS_tempmeantimeindex(GWI_IncreaseOK);
    ABS_tempmeantime = ABS_tempmeantime(GWI_IncreaseOK);
    ABS_temptimeindex = ABS_temptimeindex(GWI_IncreaseOK);
    ABS_temptime = ABS_temptime(GWI_IncreaseOK);
    ABS_tempTSmean = ABS_tempTSmean(GWI_IncreaseOK);
    ABS_tempcdpmean = ABS_tempcdpmean(GWI_IncreaseOK);
    ABS_tempBaseflowCorrection = ABS_tempBaseflowCorrection(GWI_IncreaseOK);
    disp(strcat("After review of GWI rate of change, there are ", num2str(length(ABS_tempcount)), " storms starts."))

    % This section will delete the prestorm (preoutlier groups) points if it is an outlier.
    [Outlierlocs] = isoutlier(ABS_tempBaseflowCorrection);
    outind = find(Outlierlocs==0);
    ABS_tempcount = ABS_tempcount(outind);
    ABS_tempmeantimeindex = ABS_tempmeantimeindex(outind);
    ABS_tempmeantime = ABS_tempmeantime(outind);
    ABS_temptimeindex = ABS_temptimeindex(outind);
    ABS_temptime = ABS_temptime(outind);
    ABS_tempTSmean = ABS_tempTSmean(outind);
    ABS_tempcdpmean = ABS_tempcdpmean(outind);
    ABS_tempBaseflowCorrection = ABS_tempBaseflowCorrection(outind);
    disp(strcat("After review of GWI Outliers, there are ", num2str(count), " storms starts."))

    AveBeforeStorm.count = ABS_tempcount;
    AveBeforeStorm.meantimeindex = ABS_tempmeantimeindex;
    AveBeforeStorm.meantime = ABS_tempmeantime;
    AveBeforeStorm.timeindex = ABS_temptimeindex;
    AveBeforeStorm.time = ABS_temptime;
    AveBeforeStorm.TSmean = ABS_tempTSmean;
    AveBeforeStorm.cdpmean = ABS_tempcdpmean;
    AveBeforeStorm.BaseflowCorrection = ABS_tempBaseflowCorrection;

end


%*************************************************************************



%*************************************************************************

function minsfound = Findmins
global NP
  %  ********* Begin Dynamic Warping *******************
  
    Dlocmin = islocalmin(DailyAve);
    Minimumdays = Days(Dlocmin);
    Minimumdaysval= DailyAve(Dlocmin);
    DailyavebypointMA =  zeros(1,NP);
    % looking backward 24hr average
    for k10 = 1:NP
        DailyavebypointMA(k10) = mean(Q(max(1,k10-TSperDay+1):k10));
    end

    %find 48hr minima points
    DailyavebypointMAbl =  DailyavebypointMA;
    I1 = 0;
    while min(DailyavebypointMAbl)<1000000
        I1 = I1+1;
        Curmin =find(DailyavebypointMAbl ==min(DailyavebypointMAbl),1);
        Qmins48(I1) = DailyavebypointMA(Curmin);
        Tmins48DS(I1) = DS(Curmin);
        Tmins48DN(I1) = DN(Curmin);
        DailyavebypointMAbl( max(1,Curmin-TSperDay*14) : min(NP,Curmin+TSperDay*14)) = 1000000;
    end

    [Temp,Ind] = sort(Tmins48DN)

       Qmins48 = Qmins48(Ind);
       Tmins48DN = Tmins48DN(Ind);
       Tmins48DS = Tmins48DS(Ind)

       sQmins48= Qmins48;
       sTmins48DN= Tmins48DN;
       sTmins48DS = Tmins48DS;

    %need to develope procedure to keep at least one day per month.

       TF=0;
       Threshhold2 = 7
       for k=1:10
          while sum(TF)<1
              Threshhold2 =Threshhold2-.001
              TF = isoutlier(Qmins48,'movmedian', 20, 'ThresholdFactor', Threshhold2)
          end

          figure(1) %Plot The flow and Baseflow pattern
    %       clf
    %       hold on
    %       plot(sTmins48DN,sQmins48, 'b')
    %       plot(Tmins48DN,Qmins48,'r', Tmins48DN(TF),Qmins48(TF),'gx')
    %       hold off

          Qmins48 = Qmins48(TF==0);
          Tmins48DS = Tmins48DS(TF==0);
          Tmins48DN = Tmins48DN(TF==0);

    %       figure(2) %Plot The flow and Baseflow pattern
    %       clf
    %       plot(Tmins48DN,Qmins48)
    %       
    %       figure(3) %Plot The flow and Baseflow pattern
    %       clf
    %       hold on
    %       plot(DS,Dailyavebypoint,'r')
    %       plot(DS,DailyavebypointMA,'b')
    %       plot(Tmins48DS, Qmins48,'gx')
    %       hold off
    %       pause
          TF=0;
       end

       %------- Completed Dynamic Warping -------------


end


%*************************************************************************
function [EventNums, EventStarts, EventEnds] = LoadEvents(autoload)
global DEFINED PATH
   if autoload == 0
       [filename,path] = uigetfile('*.evt')
       eventfile = fileread([path filename]);
   elseif autoload ==1
       filename = 'Storm_Events.evt';
       eventfile = fileread([PATH filename]);
   elseif autoload ==2
       filename = 'MC-SW-042_Events.evt';
       eventfile = fileread([PATH filename]);
   elseif autoload ==3
       filename = 'LM-ED-047.evt';
       eventfile = fileread([PATH filename]);
   end
    EventBlockindex = strfind(eventfile, '[Event')';
    EventNumindex = strfind(eventfile, 'Name =')';
    EventStartindex = strfind(eventfile, 'Start =')';
    EventEndindex = strfind(eventfile, 'End =')';
    
    EventNums = zeros(size(EventNumindex));
    EventStartsC = zeros(size(EventNumindex));
    EventEndsC = zeros(size(EventNumindex));
    
 
    if EventNumindex(1) < EventStartindex(1) && EventStartindex(1) < EventEndindex(1)
        for k5 = 1:length(EventNumindex)
            EventNums(k5)  = str2num(eventfile(EventNumindex(k5)+6: EventStartindex(k5)-1));
            EventStartsC(k5) = str2num(eventfile(EventStartindex(k5)+7:EventEndindex(k5)-1));
            if k5<length(EventNumindex)
                EventEndsC(k5) = str2num(eventfile(EventEndindex(k5)+5: EventBlockindex(k5+1)-2));
            else
                EventEndsC(k5) = str2num(eventfile(EventEndindex(k5)+5: length(eventfile)));
            end
        end
    
    elseif EventStartindex(1)< EventNumindex(1) && EventNumindex(1) < EventEndindex(1)
        for k5 = 1:length(EventNumindex)
            EventNums(k5)  = str2num(eventfile(EventNumindex(k5)+6: EventEndindex(k5)-1));
            EventStartsC(k5) = str2num(eventfile(EventStartindex(k5)+7:EventNumindex(k5)-1));
            if k5<length(EventNumindex)
                EventEndsC(k5) = str2num(eventfile(EventEndindex(k5)+5: EventBlockindex(k5+1)-2));
            else
                EventEndsC(k5) = str2num(eventfile(EventEndindex(k5)+5: length(eventfile)));
            end
        end
    else
        disp ('This format is not supported')
        disp ('First Number  First Start,  First End')
        disp ([ EventNumindex(1) EventStartindex(1)  EventEndindex(1)])
        return
    end
    % convertdatenumbers to Matlab
    EventStarts = datenum(EventStartsC) + datenum('1900-01-01') - 2 ; %convert to matlab epoch
    EventEnds = datenum(EventEndsC) + datenum('1900-01-01') - 2 ; %convert to matlab epoch
    datenumoffirstevent = [EventStarts(1)]   %debug
    DEFINED.Events = 1;
end







%*************************************************************************
function Cslope = Cslope_Equation(c_offset, V, Dia, depth)

   fraction = (Dia - 2 * (depth - c_offset)) / Dia;
    fraction = max(min(fraction, 1), -1);  % Clamp between -1 and 1

    arc_cos_term = 2 * acos(fraction);
    sqrt_term = sqrt(1 - fraction^2);

    R = (Dia / 8) * (1 - (2 * fraction * sqrt_term) / arc_cos_term);

    Cslope =  V / R^(2/3);
end

function c_offset = calculate_c_offset(V, Cslope, Dia, depth)
    equation = @(c_offset) offset_equation(c_offset, V, Cslope, Dia, depth);
    c_offset_initial_guess = 0.02;
    c_offset_solution = fsolve(equation, c_offset_initial_guess, optimset('Display','off'));
    c_offset = c_offset_solution;
end

function result = offset_equation(c_offset, V, Cslope, Dia, depth)
    fraction = (Dia - 2 * (depth - c_offset)) / Dia;
    fraction = max(min(fraction, 1), -1);  % Clamp between -1 and 1

    arc_cos_term = 2 * acos(fraction);
    sqrt_term = sqrt(1 - fraction^2);

    R = (Dia / 8) * (1 - (2 * fraction * sqrt_term) / arc_cos_term);

    result = R^(2/3) - (V / Cslope);
end

function [BaseTS] = F_baseTScreate(DNt,type)
    global PATTERN
 % Purpose: Generate a time series typical pattern from the seven day pattern..
    
    % Incoming Variables
        % PATTERN
            % PATTERN.num  There are at least 4 diffeent patterns that can be generated
            % PATTERN.name text string of the name of the abbreviated name of the variable
            % PATTERN.LongName  text string of the longer more clear name of the variable
            % PATTERN.Medianbyindex the particualar pattern may be one day, two
                % days or seven days.  This is the median of each index point.  In
                % general it is only used to generate the appropraite seven day
                % pattern which is stored in PATTERN.as7day.  
            % PATTERN.as7day this is the pattern converted to a seven day
                % pattern.  Since this is a seven day pattern no conversion is
                % necessary
        % DNt timedate series as numbers
    
    % Outgoing Variables
        % PATTERN as defined in Incoming variables
        % BaseTS unadjusted base time series  corresponding to DNt using type 1-4
    
    % Temporary Variables
        % Hr the hour of the day corresponding to each value in the DN the date number variable
        % MNr  the minute of the hour corresponding to each value in the DN the date number variable
        % datapntindex
        % TSatindex the seven-day pattern data point per week same as PATTERN.as7day
%******************************************************************************************************
    
    TSatindex = PATTERN(type).as7day;
    [~,~,~,Hr,MNr,~] = datevec(DNt);
    datapntindex = 1 + 288*(weekday(DNt)-1) + 12*Hr + floor(MNr/5);% Rounds down to nearest five minutes
    BaseTS = zeros(size(datapntindex));
    for k4 = 1:length(TSatindex)
        BaseTS(datapntindex==k4) = TSatindex(k4);
    end
end


function  F_DailyPatternGenerate_1day (TS, DNt)
   global PATTERN 

   % Purpose: Generate a seven day typical pattern from TS and DNt using the 1 day. This is type 4.
    
    % Incoming Variables
        % PATTERN
            % PATTERN.num  There are at least 4 diffeent patterns that can be generated
            % PATTERN.name text string of the name of the abbreviated name of the variable
            % PATTERN.LongName  text string of the longer more clear name of the variable
            % PATTERN.Medianbyindex the particualar pattern may be one day, two
                % days or seven days.  This is the median of each index point.  In
                % general it is only used to generate the appropraite seven day
                % pattern which is stored in PATTERN.as7day.  
            % PATTERN.as7day this is the pattern converted to a seven day
                % pattern.  It is generated by adding the 1 day together 7
                % times
        % TS  Time series of data, typically level , velocity or flow.
        % Typically this has already been smoothed by Savitsky-Golay
        % DNt  timedate series as numbers
    
    % Outgoing Variables
        % PATTERN as defined in Incoming variables
    
    % Temporary Variables
        % Hr the hour of the day corresponding to each value in the DN the date number variable
        % MNr  the minute of the hour corresponding to each value in the DN the date number variable
        % datapntindex
        % Medianbyindex the median value of the TS at each index defined by datapntindex
        
%******************************************************************************************************
    
    [~,~,~,Hr,MNr,~] = datevec(DNt);
    datapntindex = 1 + 288*(weekday(DNt)-1) + 12*Hr + floor(MNr/5);% Rounds down to nearest five minutes.
    datapntindex = datapntindex - 24*12*(weekday(DNt)-1);  %one day data pattern
    Medianbyindex = zeros(1,288);
    for k7 = 1:288
        Medianbyindex(k7) = median(TS(datapntindex==k7 ));%median
    end
    PATTERN(4).num = 4;
    PATTERN(4).name = 'OneDayMed';
    PATTERN(4).LongName = "1-day";
    PATTERN(4).Medianbyindex = Medianbyindex;
    PATTERN(4).as7day = [Medianbyindex Medianbyindex Medianbyindex Medianbyindex Medianbyindex Medianbyindex Medianbyindex ];
end

%*************************************************************************
function F_DailyPatternGenerate_2day (TS, DNt)
   global PATTERN


   % Purpose: Generate a seven day typical pattern from TS and DNt using a weekday pattern and a weekend pattern.  
       % The weekday starts on Sunday eveneing at 6:00 PM and has five identical days.  
       % Similarly the Weekend pattern starts on Friday at 6:00 Pm and has
       % two identical days.
    
    % Incoming Variables
        % PATTERN
            % PATTERN.num  There are at least 4 diffeent patterns that can be generated
            % PATTERN.name text string of the name of the abbreviated name of the variable
            % PATTERN.LongName  text string of the longer more clear name of the variable
            % PATTERN.Medianbyindex the particualar pattern may be one day, two
                % days or seven days.  This is the median of each index point.  In
                % general it is only used to generate the appropraite seven day
                % pattern which is stored in PATTERN.as7day.  
            % PATTERN.as7day this is the pattern converted to a seven day
                % pattern.  It is generated by combining the 2-days
        % TS  Time series of data, typically level , velocity or flow.
            % Typically this has already been smoothed by Savitsky-Golay
        % DNt  timedate series as numbers
    
    % Outgoing Variables
        % PATTERN as defined in incoming variables
    
    % Temporary Variables
        % Hr the hour of the day corresponding to each value in the DN the date number variable
        % MNr  the minute of the hour corresponding to each value in the DN the date number variable
        % datapntindex
        % Medianbyindex the median value of the TS at each index defined by datapntindex
%******************************************************************************************************
    
        DNadj = DNt-(18/24);% change start of week to 6:00 PM Sunday
        [~,~,~,Hadj,MNadj,~] = datevec(DNadj);
        Datapntofweekadj = 1 + 288*(weekday(DNadj)-1) + 12*Hadj + floor(MNadj/5);
        Wkdaywkend = zeros (size(DNt));
        Wkdaywkend(weekday(DNadj)<6)= 1;
        Wkdaywkend(weekday(DNadj)>5)= 2;
        Datapnt2day =Datapntofweekadj - 24*12*(weekday(DNadj)-1) + 288*(Wkdaywkend-1);
        Medianbyindex = zeros(1,288*2);
        for k8 = 1:288*2
            Medianbyindex(k8) = median(TS(Datapnt2day==k8 ));%median
        end
        PATTERN(3).num = 3;
        PATTERN(3).name = 'TwoDayMed';
        PATTERN(3).LongName = "2-day Weekday/Weekend";
        PATTERN(3).Medianbyindex = Medianbyindex;
        % Build the seven day pattern from the 2-day pattern
        PATTERN(3).as7day = zeros(1,2016);
            PATTERN(3).as7day(1:214) = Medianbyindex(363:576);
            PATTERN(3).as7day(215:502) = Medianbyindex(1:288);
            PATTERN(3).as7day(503:790) = Medianbyindex(1:288);
            PATTERN(3).as7day(791:1078) = Medianbyindex(1:288);
            PATTERN(3).as7day(1079:1366) = Medianbyindex(1:288);
            PATTERN(3).as7day(1367:1942) = Medianbyindex;
            PATTERN(3).as7day(1943:2016) = Medianbyindex(289:362);
end



function F_DailyPatternGenerate_diff (TS, DNt)
   global PATTERN

    % Purpose: Generate a seven day typical pattern from TS and DNt using
    % the difference method. This is type 1.
    
    % Incoming Variables
        % PATTERN
            % PATTERN.num  There are at least 4 diffeent patterns that can be generated
            % PATTERN.name text string of the name of the abbreviated name of the variable
            % PATTERN.LongName  text string of the longer more clear name of the variable
            % PATTERN.Medianbyindex the particualar pattern may be one day, two
                % days or seven days.  This is the median of each index point.  In
                % general it is only used to generate the appropraite seven day
                % pattern which is stored in PATTERN.as7day.  
            % PATTERN.as7day this is the pattern converted to a seven day
                % pattern.  Since this is a seven day pattern no conversion is
                % necessary
        % TS  Time series of data, typically level , velocity or flow.
        % Typically this has already been smoothed by Savitsky-Golay
        % DNt  timedate series as numbers
    
    % Outgoing Variables
        % PATTERN as defined in Incoming variables
    
    % Temporary Variables
        % Hr the hour of the day corresponding to each value in the DN the date number variable
        % MNr  the minute of the hour corresponding to each value in the DN the date number variable
        % datapntindex
        
        % datapntofweekdiff used for data points of differencing
        % Qdiff the difference of the time series. one point less than TS
        % Mediandiffbyindex median difference for each point in the week.
        % Medianfromdiffbyindex this is the daily pattern made by crawling along the individual median differences
        % Scale used to reset the slope of the difference daily pattern
        % Medianbyindex the median value of the TS at each index defined by datapntindex
        % linfit should equal zero if the differences are not climbing or falling
        % Slopdif just a text string of th eslope to be displayed if not level
        
%******************************************************************************************************
    [~,~,~,Hr,MNr,~] = datevec(DNt);
    datapntindex = 1 + 288*(weekday(DNt)-1) + 12*Hr + floor(MNr/5);% Rounds down to nearest five minutes.

        datapntofweekdiff = datapntindex(1:length(datapntindex)-1); % use for data points of differencing
        Qdiff = diff(TS); % for the differencing method
        Mediandiffbyindex = zeros(1,2016);
        for k3 =1:2016
            Mediandiffbyindex(k3) = median(Qdiff(datapntofweekdiff==k3 ));%median diff
        end
        Medianfromdiffbyindex = zeros(1,2016);
        Medianfromdiffbyindex(1) = median(TS(datapntindex==1));
        for k4 = 2:2016
            Medianfromdiffbyindex(k4) = Medianfromdiffbyindex(k4-1)+ Mediandiffbyindex(k4-1);
        end
        
        %Set daily pattern from diff to level
        Scale = (Medianfromdiffbyindex(1) -  (Medianfromdiffbyindex(2016)+ Mediandiffbyindex(2016)))/2016;
        Medianbyindex  = Medianfromdiffbyindex +(0:2015)*Scale;
        
        %Check daily pattern is level
        linfit = polyfit(1:2016,Medianbyindex,1);
        Slopdif = num2str(round(linfit(1)),4);
        disp(" ")
        disp(strcat("Slope of line through Diff Points = ", Slopdif))
        if Medianfromdiffbyindex(1) ~=   Medianbyindex(2016)+Mediandiffbyindex(2016)+Scale
            disp('Warning unbalanced on 7-day')
            linfit = polyfit(1:2016,Medianbyindex,1);
            disp(' ')
            disp('Use keyboard function to investigate why it is not balanced.')
            %keyboard
        else
            disp("Balanced")
        end

        %reassign weekly diff pattern to a baseflow
        PATTERN(1).num = 1;
        PATTERN(1).name = 'SevDayDif'; 
        PATTERN(1).LongName = "7-day by Differencing";
        PATTERN(1).Medianbyindex = Medianbyindex;
        PATTERN(1).as7day = Medianbyindex; 
end


function F_DailyPatternGenerate_Median (TS, DNt)
global PATTERN
    % Purpose: Generate a seven day typical pattern from TS and DNt using the median method. This is type 2.
    
    % Incoming Variables
        % PATTERN
            % PATTERN.num  There are at least 4 diffeent patterns that can be generated
            % PATTERN.name text string of the name of the abbreviated name of the variable
            % PATTERN.LongName  text string of the longer more clear name of the variable
            % PATTERN.Medianbyindex the particualar pattern may be one day, two
                % days or seven days.  This is the median of each index point.  In
                % general it is only used to generate the appropraite seven day
                % pattern which is stored in PATTERN.as7day.  
            % PATTERN.as7day this is the pattern converted to a seven day
                % pattern.  Since this is a seven day pattern no conversion is
                % necessary
        % TS  Time series of data, typically level , velocity or flow.
        % Typically this has already been smoothed by Savitsky-Golay
        % DN  timedate series as numbers
    
    % Outgoing Variables
        % PATTERN as defined in Incoming variables
    
    % Temporary Variables
        % Hr the hour of the day corresponding to each value in the DN the date number variable
        % MNr  the minute of the hour corresponding to each value in the DN the date number variable
        % datapntindex
        % Medianbyindex the median value of the TS at each index defined by datapntindex
        
%******************************************************************************************************

    [~,~,~,Hr,MNr,~] = datevec(DNt);
    datapntindex = 1 + 288*(weekday(DNt)-1) + 12*Hr + floor(MNr/5);% Rounds down to nearest five minutes
    Medianbyindex = zeros(1,2016);
    for k3 =1:2016
        Medianbyindex(k3) = median(TS(datapntindex==k3));%median at each point in the week
    end
    PATTERN(2).num = 2;
    PATTERN(2).name = 'SevDayMed';
    PATTERN(2).LongName = "7-day by Median";
    PATTERN(2).Medianbyindex = Medianbyindex;
    PATTERN(2).as7day = Medianbyindex;

end















